Cholesky factorization as structural factorization (cholsvar.prg)

The following program illustrates the two forms of imposing the identifying restrictions in structural decompositions. For illustration purposes and to check that the restrictions are correctly imposed, we impose restrictions that replicate the Cholesky factorization. (If all you want is the standard impulse response based on the Cholesky factorization, there is no need to use the structural decomposition.) The program checks whether the factorization matrix from the three methods match by computing the L-infinity norm (maximum of absolute column sums) of the matrix difference. Theoretically, the difference should be zero but because of floating errors the numerical difference won't be exactly zero.

 
' replicate cholesky factorization using SVAR
' 11/5/99 h
' last checked 3/7/2007
 
' include subroutine
include sub_rmaxdiff.prg
 
'change path to program path
%path = @runpath
cd %path
 
' create workfile
wfcreate cholsvar q 1948:1 1979:3
 
' fetch data from database
fetch(d=data_svar) rgnp rinv m1
 
' take log of each series
series lgnp = log(rgnp)
series linv = log(rinv)
series lm1 = log(m1)
 
' estimate unrestricted VAR
var var1.ls 1 4 lgnp linv lm1 @ c
 
'-------------------------------------------------------------------
' method 1: short-run restrictions in text form
'-------------------------------------------------------------------
 
var1.cleartext(svar)
var1.append(svar) @e1 = c(1)*@u1
var1.append(svar) @e2 = -c(2)*@e1 + c(3)*@u2
var1.append(svar) @e3 = -c(4)*@e1 - c(5)*@e2 + c(6)*@u3
 
freeze(tab1) var1.svar(rtype=text,conv=1e-5)
'show tab1
 
' store estimated A and B matrices from method 1
matrix mata1 = var1.@svaramat
matrix matb1 = var1.@svarbmat
 
' compute factorization matrix
matrix fact1 = @inverse(mata1)*matb1
 
'-------------------------------------------------------------------
' method 2: short-run restrictions in pattern matrix
'-------------------------------------------------------------------
 
' create pattern matrices
matrix(3,3) pata
' fill matrix in row major order
pata.fill(by=r) 1,0,0, na,1,0, na,na,1
 
matrix(3,3) patb
patb.fill(by=r) na,0,0, 0,na,0, 0,0,na
 
freeze(tab2) var1.svar(rtype=patsr,conv=1e-5,namea=pata,nameb=patb)
'show tab2
 
' store estimated A and B matrices from method 2
matrix mata2 = var1.@svaramat
matrix matb2 = var1.@svarbmat
 
' compute factorization matrix
matrix fact2 = @inverse(mata2)*matb2
 
'-------------------------------------------------------------------
' method 3: built-in cholesky
'-------------------------------------------------------------------
 
' do impulse response with default Cholesky ordering
var1.impulse(10)
 
' store cholesky factor
matrix fact3 = var1.@impfact
 
'-------------------------------------------------------------------
' check whether results match
'-------------------------------------------------------------------
 
' compute difference in factorization matrices
scalar diff13
call sub_rmaxdiff(fact1, fact3, diff13)
scalar diff23
call sub_rmaxdiff(fact2, fact3, diff23)
 
' display L-infinity norm of matrix difference in table
table(2,2) check
setcolwidth(check,1,20)
check(1,1) = "method1 - method3:"
check(1,2) = diff13
check(2,1) = "method2 - method3:"
check(2,2) = diff23
 
show check
 

^return

Blanchard-Quah (1989) long-run restriction (blanquah1.prg)

The following program replicates the long-run restriction model by Blanchard-Quah (1989) (the data are not the same as those in the original so that the results do not match) .The DY (output growth) and U (unemployment) series are already demeaned following the procedure by Blanchard-Quah (1989). The program estimates the factorization matrix by three different methods which should all yield the same result.

' Blanchard-Quah long-run restriction (11/5/99)
' verify estimates using three methods
' last checked 3/7/2007
 
' include subroutine
include sub_rmaxdiff.prg
 
'change path to program path
%path = @runpath
cd %path
 
' create workfile
wfcreate blanquah q 1948:1 1987:4
 
' fetch data from database (already demeaned)
fetch(d=data_svar) dy u
 
' estimate VAR (no constant)
var var1.ls 1 8 dy u @ 
 
'-------------------------------------------------------------------
' method 1: text form long-run restrictions
'-------------------------------------------------------------------
 
var1.cleartext(svar)
var1.append(svar) @lr1(@u1)=0
 
freeze(tab1) var1.svar(rtype=text,conv=1e-4)
show tab1
 
' store estimated A and B matrices from method 1
matrix mata1 = var1.@svaramat    ' A should be identity
matrix matb1 = var1.@svarbmat
 
'-------------------------------------------------------------------
' method 2: long-run restrictions in short-run form 
' *not recommended; only for checking*
'-------------------------------------------------------------------
 
' get unit (cumulated) long-run response
var1.impulse(imp=u)
matrix clr = var1.@lrrsp
 
var1.cleartext(svar)
var1.append(svar) @e1 = c(1)*@u1 + c(2)*@u2
var1.append(svar) @e2 = -clr(1,1)*c(1)/clr(1,2)*@u1 + c(4)*@u2
 
freeze(tab2) var1.svar(rtype=text,conv=1e-5)
show tab2
 
' store estimated A and B matrices from method 2
matrix mata2 = var1.@svaramat    ' A should be identity
matrix matb2 = var1.@svarbmat
 
'-------------------------------------------------------------------
' method 3: moment matching (used in Blanchard-Quah (1989) paper) 
' *only works for just-identified models*
'-------------------------------------------------------------------
 
' get residual covariance 
sym rcov = var1.@residcov
' and do cholesky factorization
matrix chol = @cholesky(rcov)
 
' factorize unit long-run response
matrix p = clr * chol
 
' factorized B matrix 
matrix(2,2) q
' impose sign normalization by taking positive root
q(1,1) = @sqrt( 1/(1+p(1,1)*p(1,1)/p(1,2)/p(1,2)) )
q(2,1) = -p(1,1)*q(1,1)/p(1,2)
q(1,2) = -q(2,1)
q(2,2) = q(1,1)
 
' get implied B matrix
matrix matb3 = chol * q
 
'-------------------------------------------------------------------
' check whether results match
'-------------------------------------------------------------------
 
' difference in A matrices
matrix eye2 = @identity(2)       ' truth
scalar adiff1
call sub_rmaxdiff(mata1, eye2, adiff1)
scalar adiff2
call sub_rmaxdiff(mata2, eye2, adiff2)
 
' difference in B matrices
scalar bdiff13
call sub_rmaxdiff(matb1, matb3, bdiff13)
scalar bdiff23
call sub_rmaxdiff(matb2, matb3, bdiff23)
 
' display L-infinity norm of matrix difference in table
table(2,3) check
setcolwidth(check,1,10)
setcolwidth(check,2,15)
setcolwidth(check,3,15)
 
check(1,2) = "A - eye2"
check(1,3) = "B - method3"
check(2,1) = "method1:"
check(3,1) = "method2:"
 
check(2,2) = adiff1
check(3,2) = adiff2
check(2,3) = bdiff13
check(3,3) = bdiff23
 
show check
 

^return

Short-run restriction impulse responses (Sims 1986) (sims1.prg)

The following program replicates the short-run restriction exercises in Sims (1986) (the data are not the same as those in the original so that the results do not match). Note that while Sims assumes a diagonal covariance matrix for the structural innovations, EViews assumes an identity covariance matrix. For this reason, we need to estimate the standard deviation of the structural shocks as elements of the B matrix.

 
' replicate Sims (1986) with similar data 
' 1/12/2000 h
' last checked 3/7/2007
 
'change path to program path
%path = @runpath
cd %path
 
' create workfile                                                                       
wfcreate sims1 q 1948:1 1979:3
 
' fetch data from database
fetch(d=data_svar) rgnp rinv pidgnp m1 tb3sec unemp
 
' transform series 
series lgnp = log(rgnp)
series linv = log(rinv)
series lprc = log(pidgnp)
series lm1 = log(m1)
 
'----------------------------------------------------------------------
' replicate Chart 1 (p.11) *no Bayesian priors imposed*
'----------------------------------------------------------------------
 
var var1.ls 1 4 lgnp linv lprc lm1 unemp tb3sec @ c
freeze(chart1) var1.impulse(32,m,imp=chol,se=a)
show chart1
 
' plot only some of money innovation responses
freeze(chart1a) var1.impulse(32,m,imp=chol,se=a) lgnp lprc tb3sec @ lm1 
 
'----------------------------------------------------------------------
' replicate Chart 2 (p.13)
'----------------------------------------------------------------------
 
' reorder VAR to confirm to Sims' (1986) ordering
var var1.ls 1 4 tb3sec lm1 lgnp lprc unemp linv @ c
 
' 1st identification restrictions (p.12 (12)-(17))
coef(6) b
var1.append(svar) @e1 = c(1)*@e2 + b(1)*@u1
var1.append(svar) @e2 = c(2)*@e3 + c(3)*@e4 + c(4)*@e1 + b(2)*@u2
var1.append(svar) @e3 = c(5)*@e1 + c(6)*@e6 + b(4)*@u4
var1.append(svar) @e4 = c(7)*@e1 + c(8)*@e3 + c(9)*@e6 + b(5)*@u5
var1.append(svar) @e5 = c(10)*@e1 + c(11)*@e3 + c(12)*@e6 + c(13)*@e4 + b(6)*@u6
var1.append(svar) @e6 = b(3)*@u3
 
' estimate structural factorization
freeze(svar1) var1.svar(rtype=text,conv=1e-5)
show svar1
 
' replicate Chart 2 (p.13)
freeze(chart2) var1.impulse(32,m,imp=struct,se=a) lgnp linv lprc lm1 unemp tb3sec @ 1 2 4 5 6 3
show chart2
 
' plot only some of money supply responses
freeze(chart2a) var1.impulse(32,m,imp=struct,se=a) lgnp lprc tb3sec @ 1 
 
'----------------------------------------------------------------------
' replicate Chart 3 (p.13)
'----------------------------------------------------------------------
 
' clear previous restrictions
var1.cleartext(svar)
 
' 2nd identification restrictions (p.14 (18)-(23))
var1.append(svar) @e1 = c(1)*@e2 + b(1)*@u1
var1.append(svar) @e2 = c(2)*@e3 + c(3)*@e6 + c(4)*@e4 + c(5)*@e1 + b(2)*@u2
var1.append(svar) @e3 = c(6)*@e1 + c(7)*@e6 + b(3)*@u3
var1.append(svar) @e4 = c(8)*@e1 + c(9)*@e3 + c(10)*@e2 + b(5)*@u5
var1.append(svar) @e5 = c(11)*@e1 + c(12)*@e3 + c(13)*@e4 + c(14)*@e6 + b(6)*@u6
var1.append(svar) @e6 = b(4)*@u4
 
' estimate structural factorization
' default starting value does not converge
rndseed 12
freeze(svar2) var1.svar(rtype=text,conv=1e-5,f0=u)
show svar2
 
' replicate Chart 3 (p.13)
freeze(chart3) var1.impulse(32,m,imp=struct,se=a) lgnp linv lprc lm1 unemp tb3sec @ 1 2 3 5 6 4
show chart3
 
' plot only some of money supply responses
freeze(chart3a) var1.impulse(32,m,imp=struct,se=a) lgnp lprc tb3sec @ 1
 
 

^return

Long-run restriction impulse responses (Blanchard-Quah 1989) (blanquah2.prg)

The following program replicates the long-run restriction model by Blanchard-Quah (1989) (the data are not the same as those in the original so that the results do not match). The DY (output growth) and U (unemployment) series are already demeaned following the procedure in Blanchard-Quah (1989). To plot the accumulated response and the level response in one graph, we first store the impulse responses in separate matrices, extract and place the relevant columns into a third matrix, and use the line view of that matrix.

 
' replicate impulse response graph
' from long-run restriction (Blanchard-Quah 1989)
' 11/12/99 h
' checked 3/7/2007
 
' include subroutine
include sub_rmaxdiff.prg
 
'change path to program path
%path = @runpath
cd %path
 
' create workfile
wfcreate blanquah q 1948:1 1987:4
 
' fetch data from database (already demeaned)
fetch(d=data_svar) dy u
 
' change to percentage
dy = 100*dy
 
' estimate VAR (no constant)
var var1.ls 1 8 dy u @ 
 
' impose long-run restrictions
' @u1 = aggregate demand shock
' @u2 = aggregate supply shock
var1.cleartext(svar)
var1.append(svar) @lr1(@u1)=0
 
' estimate factorization
freeze(tab1) var1.svar(rtype=text,conv=1e-5)
 
' set response horizon
!hrz = 40
 
' replicate Figures 3-4 (p.662)
' accumulate to get level response of output
freeze(fig3) var1.impulse(!hrz,imp=struct,se=a,a,matbys=rsp_acc) dy @ 1
freeze(fig4) var1.impulse(!hrz,imp=struct,se=a,a) dy @ 2
freeze(fig34) fig3 fig4
show fig34
 
' replicate Figures 5-6
freeze(fig5) var1.impulse(!hrz,imp=struct,se=a,matbys=rsp_lvl) u @ 1
freeze(fig6) var1.impulse(!hrz,imp=struct,se=a) u @ 2
freeze(fig56) fig5 fig6
show fig56
 
' replicate Figure 1
matrix(!hrz,2) mtmp
' get accumulated output response
vector v = @columnextract(rsp_acc,1)
colplace(mtmp,v,1)
' get level unemployment response
v = @columnextract(rsp_lvl,2)
colplace(mtmp,v,2)
 
freeze(fig1) mtmp.line
fig1.draw(line,left) 0
fig1.elem(1) legend(Output response to demand)
fig1.elem(2) legend(Unemployment response to demand)
show fig1
 
' replicate Figure 2
' get accumulated output response
vector v = @columnextract(rsp_acc,3)
colplace(mtmp,v,1)
' get level unemployment response
v = @columnextract(rsp_lvl,4)
colplace(mtmp,v,2)
 
freeze(fig2) mtmp.line
fig2.draw(line,left) 0
fig2.elem(1) legend(Output response to supply)
fig2.elem(2) legend(Unemployment response to supply)
show fig2
 
 

^return

Bootstrap impulse response confidence bounds for structural VARs with long-run restrictions (bqboot.prg)

EViews does not currently provide standard error bounds for impulse responses from structural impulses identified by long-run restrictions. The following program illustrates how to obtain bootstrap impulse response confidence bounds for the Blanchard-Quah (1989) model. The bootstrap is performed as follows:

1. Resample from the VAR residuals with replacement.
2. Obtain bootstrap data using the estimated VAR coefficients and the bootstrap residuals.
3. Estimate the VAR and structural factorization using the bootstrap data from step 2.
4. Store the impulse responses using estimates from step 3.
5. Repeat steps 1-4 many times.

Step 2 is implemented by specifying the bootstrap residuals as (level shift) add factors in model solve. Note also that estimation of the structural factorization at each bootstrap step may fail to converge; in such cases the code only stores results if convergence was achieved. The code below performs 100 bootstrap replications; you should probably increase this number to get more reliable estimates but be aware that this make take a very long time to execute.

 
' bootstrap structural VAR impulse responses (11/28/2000 )
' from long-run restriction (Blanchard-Quah 1989)
' last checked 3/7/2007
 
include sub_bootirf.prg
 
'change path to program path
%path = @runpath
cd %path
 
' create workfile
wfcreate blanquah q 1948:1 1987:4
 
' fetch data from database (already demeaned)
fetch(d=data_svar) dy u
 
' change to percentage
dy = 100*dy
 
' estimate VAR (no constant)
var var1.ls 1 8 dy u @ 
 
' impose long-run restrictions
' @u1 = aggregate demand shock
' @u2 = aggregate supply shock
var1.cleartext(svar)
var1.append(svar) @lr1(@u1)=0
 
' estimate factorization
var1.svar(rtype=text,conv=1e-5)
 
' set response horizon
!hrz = 40
 
' store level and accumulated responses
do var1.impulse(!hrz,imp=struct,matbys=rsp_lvl) dy u @ 1
do var1.impulse(!hrz,imp=struct,a,matbys=rsp_acc) dy u @ 1
 
' get residuals to bootstrap
var1.makeresid(n=gres) res_dy res_u
 
' make model to generate bootstrap data
var1.makemodel(mod1)
' assign add factors for bootstrap residuals
mod1.addassign(i) @all
 
'set random number generator
rndseed(type=mt) 1234567
!reps = 100
'declare storage matrices
for !i=1 to @columns(rsp_lvl)
   matrix(!hrz,!reps) brsp_lvl{!i}
   matrix(!hrz,!reps) brsp_acc{!i}
next
'temporary working vector
vector vtmp
'count number of cases that converged
scalar cnt = 0
 
'bootstrap loop
for !i=1 to !reps
   'bootstrap residuals
   smpl @all if res_dy<>na
   gres.resample dy_a u_a
   'generate bootstrap data
   mod1.solve
   'estimate VAR (no constant) with bootstrap data
   smpl @all
   var var2.ls 1 8 dy_0 u_0 @ 
   'impose long-run restrictions
   var2.cleartext(svar)
   var2.append(svar) @lr1(@u1)=0
   'estimate factorization
   var2.svar(rtype=text,conv=1e-5,maxiter=100,nostop)
   'store results only if converged
   if (var2.@svarcvgtype=0) then
          'increment counter
          cnt = cnt + 1
          'get level and accumulated bootstrap responses
          do var2.impulse(!hrz,imp=struct,matbys=tmp_lvl) dy_0 u_0 @ 1
          do var2.impulse(!hrz,imp=struct,a,matbys=tmp_acc) dy_0 u_0 @ 1
          'and store in appropriate position
          for !j=1 to @columns(rsp_lvl)
                  vtmp = @columnextract(tmp_lvl,!j)
                  colplace(brsp_lvl{!j}, vtmp, cnt)
                  vtmp = @columnextract(tmp_acc,!j)
                  colplace(brsp_acc{!j}, vtmp, cnt)
          next
          delete tmp_lvl tmp_acc
   endif
next
 
'reshape results matrix if there were non-converged cases
if (cnt <> !reps) then
   for !i=1 to @columns(rsp_lvl)
          brsp_lvl{!i} = @subextract(brsp_lvl{!i}, 1, 1, !hrz, cnt)
          brsp_acc{!i} = @subextract(brsp_acc{!i}, 1, 1, !hrz, cnt)
   next
   table warning
   %msg = @str(!reps-cnt) + " case(s) out of " + @str(!reps) + " replications failed to converge"
   warning(1,1) = %msg
   show warning
endif
 
'find bootstrap percentiles for each horizon
matrix(!hrz,@columns(rsp_lvl)) lvl_upp
matrix(!hrz,@columns(rsp_lvl)) lvl_low
matrix(!hrz,@columns(rsp_lvl)) acc_upp
matrix(!hrz,@columns(rsp_lvl)) acc_low
 
rowvector rtmp
for !i=1 to !hrz
   for !j=1 to @columns(rsp_lvl)
          rtmp = @rowextract(brsp_lvl{!j},!i)
          lvl_upp(!i,!j) = @quantile(rtmp,0.833)
          lvl_low(!i,!j) = @quantile(rtmp,0.167)
          rtmp = @rowextract(brsp_acc{!j},!i)
          acc_upp(!i,!j) = @quantile(rtmp,0.833)
          acc_low(!i,!j) = @quantile(rtmp,0.167)
   next
next
 
'plot impulse responses with bootstrap error bounds
%gname = "fig_lvl"
call sub_bootirf(rsp_lvl, lvl_upp, lvl_low, %gname)
' add zero line to all graphs
{%gname}.draw(line,left,rgb(180,180,180)) 0
' add title to graph
{%gname}.addtext(0.1,-0.2) Bootstrap 66% Confidence Bounds (Levels)
' realign
{%gname}.align(2,1,0.5)
show {%gname}
 
%gname = "fig_acc"
call sub_bootirf(rsp_acc, acc_upp, acc_low, %gname)
' add zero line to all graphs
{%gname}.draw(line,left,rgb(180,180,180)) 0
' add title to graph
{%gname}.addtext(0.1,-0.2) Bootstrap 66% Confidence Bounds (Accumulated Responses)
' realign
{%gname}.align(2,1,0.5)
show {%gname}
 
 

^return